
/*
 *	round_geom.c (based on geom_3d_v6.c)
 *
 *	designed for trumpet modeling with a circular bell
 *
 *	add sites L and R to represent the left and right lip sites
 *
 *	set up a nonuniform grid in 3d for use in recorder8 and later programs
 *
 *	usage: round_geom <in_file> <trumpet_geometry.dat> <output_directory> 
 *
 *	where in_file contains data used to define the nonuniform grid
 *
 *	trumpet_geometry.dat defines the instrument and is similar to (but with a few
 *	more options as) the format used by recorder_v8 and grid_3d_new_v2
 *
 *	solid points/regions with circular symmetry about an axis along x are indicated
 *		by C and the center of symmetry for that y-z row is denoted by O (capital "oh")
 *
 *	output_directory contains all the output files including the vtk files that contain
 *		the overall geometry
 *
 *	the output files are 
 *		x_grid.dat -- contains grid information for the x direction in the form
 *				0	x(0)	dx(0)
 *				1	x(1)	dx(1)
 *				etc
 *				where x(i) is the location along x of grid point i
 *				and dx(i) is the spacing between grid points i and i+1
 *		y_grid.dat -- same info for the grid along y
 *		z_grid.dat -- same info for the grid along z
 *		geometry_solid.dat -- locations of all solid points (recorder and boundary)
 *			in terms of i,j grid locations
 *		geometry_initialize.dat -- locations of the initialization region
 *			in terms of i,j grid locations
 *		geometry_sound.dat -- locations of the points where the sound is to be recorded
 *			in terms of i,j grid locations
 *
 *		geometry_instrument.dat -- locations of all solid points (recorder only)
 *			in terms of i,j grid locations
 *
 *		place recorder_geometry.dat into the output_directory for further ref
 *
 *		produce eps files showing the precise geometry of the simulation - one set of files
 *		low_res_geom_k.eps is a low res plot and high_res_geom_k.eps shows the detail in the region
 *		where the grid spacing is smallest -- the index k is the layer number as used
 *		in recorder_geometry.dat
 *
 *		the file geometry.vtk is generated for viewing with paraview
 *		and file bare_geometry.vtk is generated for viewing with paraview
 *			(this file assumes a uniform grid)
 *
 *		lip_sites.dat -- shows the layer containing the center x-y (i-j) plane of the lips
 *			indicating the location of the lip centers, the pivot points (point A
 *			in the Adachi model of lips), and ghost sites 
 *			Some important lip dimensions (in units of a grid spacing) are given
 *			at the start of this file
 *			this info is also contained in the log_file
 *
 *	-------------------------
 *
 *	the way the grid spacing varies for a given direction (i.e., the x direction) is 
 *	defined by several quantities
 *
 *	x_1 <	x_2 <	x_3 <	x_4
 *
 *	x_1: the lowest (starting) index value and dx is large at this end
 *	x_2: dx falls to its minimum value dx_min here and is constant until
 *	x_3: dx is constant between x_2 and x_3. dx grows again for larger values of x
 *	x_4: the largest (final) index value and dx is large at this end
 *
 *	the corresponding array indices are i_1, i_2, i_3, and i_4
 *
 *	add similar code for grid along y direction
 *
 *	this program reads in values for x_1, x_2, .... y_1, y_2 ... from <in_file>
 *
 *	format for in_file is
 *	x_1	x_2	x_3	x_4
 *	dx_min	dx_max	ax
 *	y_1	y_2	y_3	y_4
 *	dy_min	dy_max	ay
 *	z_1	z_2	z_3	z_4
 *	dz_min	dz_max	az
 *
 *	ax, ay, and az are factors that controls how fast dx, dy, and dz vary (higher values give a faster
 *		variation; try a = 20 first)
 *
 *	as noted above, files x_grid.dat and y_grid.dat contain 
 *		i1	x1	dx1	
 *		i2	x2	dx2	
 *		etc
 *	where i1 is the grid index
 *
 *	The files geometry_solid.dat, geometry_initialize.dat, and geometry_record.dat contain
 *	the grid coordinates i,j,k of the solid regions (boundary and recorder), the region where u
 *	is initialized, and the points where the sound is to be recorded
 *
 *	geometry_initialize.dat also contains the scaled values of the fluid speed at each point
 *		on a scale where max = 1.0
 *	
 *	The files x_grid.dat, y_grid.dat, and z_grid.dat  will be used by recorder8 to set up the computational
 *	grid, and the files geometry_solid.dat, geometry_initialize.dat, and geometry_record.dat
 *	will be used to initialize the geometry of the problem
 *
 *	also put in_file into the output directory along with some other files that will not be used
 *		but are useful for record keeping
 *
 *	the file u_init.dat contains info for initializing u in the channel
 */

#include <stdio.h>
#include <math.h>
#include <strings.h>
#include <string.h>
#include <ctype.h>
#include <stdlib.h>
#include <sys/types.h>
#include <sys/stat.h>

//#define N_X_MAX	2510
//#define N_Y_MAX	410
//#define N_Z_MAX	410

#define N_X_MAX	6010
#define N_Y_MAX	1010
#define N_Z_MAX	1010

#define YES			0
#define NO			1

#define ON			2
#define OFF			3

#define	TMP_OPEN		9
#define	OPEN			10
#define LOWER_X_BOUNDARY	11
#define UPPER_X_BOUNDARY	12
#define LOWER_Y_BOUNDARY	13
#define UPPER_Y_BOUNDARY	14
#define LOWER_Z_BOUNDARY	15
#define UPPER_Z_BOUNDARY	16
#define SOLID			17
#define INIT_POINT		18
#define MARKED			19

#define NOT_STARTED		0
#define STARTED			1

#define START_FRESH		0
#define RE_START		1

#define U			0
#define V			1
#define W			2
#define RHO			3

/* used by vector field graphing routine	*/
#define RIGHT			0
#define LEFT			1
#define UP			2
#define DOWN			3
#define UPPER_RIGHT		4
#define UPPER_LEFT		5
#define LOWER_RIGHT		6
#define LOWER_LEFT		7

/* used to classify points on the boundary of the recorder	*/
#define XUP			20
#define XDOWN			21
#define YUP			22
#define YDOWN			23
#define ZUP			24
#define ZDOWN			25
#define CORNER			26

// for lip-related points
#define	PIVOT			30
#define ORIGIN			31
#define CURRENT_CENTER		32
#define LEFT_LIP		33
#define RIGHT_LIP		34
#define GHOST			35

// for regions with circular symmetry
#define CIRCULAR		40
#define CENTER			41

#define CURRENT			0
#define OLD			1
#define NEW			2

#define	LEADING_EDGE		0
#define	TRAILING_EDGE		1
#define MIDDLE			2

#define FIRST_CIRCLE		0
#define LAST_CIRCLE		1

#define	SOUND			0
#define	VELOCITY		1

char ***allocate_array_of_chars();

double exp();
double stretch_factor();

double y_grid[N_Y_MAX];
double dy[N_Y_MAX];
double x_grid[N_X_MAX];
double dx[N_X_MAX];
double z_grid[N_Z_MAX];
double dz[N_Z_MAX];
double ax1,ax2,dx_min,dx_max;
double ay1,ay2,dy_min,dy_max;
double az1,az2,dz_min,dz_max;
double x1,x2,x3,x4;
double y2,y3,y4,y_1;
double z1,z2,z3,z4;
FILE *fp_grid_input;
FILE *fp_x_grid_dat,*fp_y_grid_dat,*fp_z_grid_dat;
FILE *fp_recorder;
double y_max,x_max,z_max;
int i_min,i_max;
int j_min,j_max;
int k_min,k_max;
int i_max_bare,j_max_bare,k_max_bare;
int i_init_min,i_init_max;
int j_init_min,j_init_max;
int k_init_min,k_init_max;
int i_init_min_bare,i_init_max_bare;
int j_init_min_bare,j_init_max_bare;
int k_init_min_bare,k_init_max_bare;

// this needs work
double init_val_bare[N_Y_MAX];
double init_val[N_Y_MAX];
int i_init_bare;
int j_init_bare_start,j_init_bare_end;

// int geometry[N_X_MAX][N_Y_MAX][N_Z_MAX]; /* describes the domain of the calc on NONUNIFORM grid	*/
// char geometry[N_X_MAX][N_Y_MAX][N_Z_MAX]; /* describes the domain of the calc on NONUNIFORM grid	*/
char ***geometry;		/* describes the domain of the calc on NONUNIFORM grid	*/
                                /* the geometry (including the shape of the recorder)   */
                                /* is specified in recorder_geometry.dat                */
                                /* most geometry-related parameters are set there       */
                                /* OPEN = open space away from walls    */
                                /* SOLID = body of recorder             */
				// LOWER_X_BOUNDARY        boundaies of simulation region
				// UPPER_X_BOUNDARY        
				// LOWER_Y_BOUNDARY       
				// UPPER_Y_BOUNDARY      
				// INIT_POINT              points where u is set

// int bare_geometry[N_X_MAX][N_Y_MAX][N_Z_MAX]; 
// char bare_geometry[N_X_MAX][N_Y_MAX][N_Z_MAX]; 
char ***bare_geometry; 		/* describes the domain of the calc on UNIFORM grid	*/

//int i_record_1,i_record_2,i_record_3;
//int j_record_1,j_record_2,j_record_3;
//int k_record_1,k_record_2,k_record_3;

//int i_record_1_bare,i_record_2_bare,i_record_3_bare;
//int j_record_1_bare,j_record_2_bare,j_record_3_bare;
//int k_record_1_bare,k_record_2_bare,k_record_3_bare;

int n_record_points;
int i_record[30],j_record[30],k_record[30];
int n_record_flag[30];	// = SOUND if only record sound, = VELOCITY if record sound and velocity

int i_record_bare[30],j_record_bare[30],k_record_bare[30];

int i_recorder_corner,j_recorder_corner,k_recorder_corner;
int i_recorder_bare,j_recorder_bare,k_recorder_bare;// bare coordinates of the lower left corner of recorder
		// WARNING - in previous versions (with uniform grids) these were the coordinates
 		// of the lower left corner
int i_init,j_init_start,j_init_end,k_init_start,k_init_end;
double x_recorder_bare,y_recorder_bare,z_recorder_bare;

double dx_bare;	// original "bare" dx for defining the recorder geometry
FILE *fp_log;

char output_directory[200],recorder_geometry[200],input_file[200];
char vtk_directory[200];

FILE *fp_diag;

double v_scale_factor_bare[N_Z_MAX];	/* scaling factors for v_x in each layer that contains an init point in bare coordinates */
double v_scale_factor[N_Z_MAX];		/* scaling factors for v_x in each layer that contains an init point in final coordinates*/

FILE *fp_lip_info;
char lip_info[200];

int
main(int argc,char *argv[])
{
	int i,j;
	char mes[500];

	fp_diag = fopen("diag_info","w");

	if(argc < 4) {
		fprintf(stderr,"usage: geometry_setup <in_file> <recorder_geometry.dat> <output_directory>\n");
		fprintf(stderr,"\tformat for input file is\n");
		fprintf(stderr,"\tx_1\tx_2\tx_3\tx_4\n");
		fprintf(stderr,"\tdx_min\tdx_max\tax\n");
		fprintf(stderr,"\ty_1\ty_2\ty_3\ty_4\n");
		fprintf(stderr,"\tdy_min\tdy_max\tay\n");
		fprintf(stderr,"\tz_1\tz_2\tz_3\tz_4\n");
		fprintf(stderr,"\tdz_min\tdz_max\taz\n");
		fprintf(stderr,"\n");
 		fprintf(stderr,"\ta(x,y,z) (= factor that controls how fast delta varies (try a = 20 first)\n");
 		fprintf(stderr,"\\ttlarger values of ax give a faster variation of dx, etc.\n");
		exit(0);
	}

	fp_log = fopen("log_file","w");
	strcpy(recorder_geometry,argv[2]);
	strcpy(output_directory,argv[3]);
	mkdir(output_directory,0777);
	sprintf(mes,"cp %s %s/",argv[1],output_directory);
	system(mes);
	sprintf(mes,"cp %s %s/",recorder_geometry,output_directory);
	system(mes);
	strcpy(vtk_directory,output_directory);
	

fprintf(stderr,"setting up grid\n");
	set_up_grid(argv[1]);
fprintf(stderr,"reading in and placing recorder\n");
	place_recorder();
	save_geometry();
//	display_geometry(1);
//	make_bare_vtk_file();
	make_vtk_file1();
//	make_vtk_file2();

	check_init_region();

	fclose(fp_diag);

	return(0);
}

int
set_up_grid(char *in_file)
{
	int i,j,k;
	char mes[500];
	FILE *fp_dx,*fp_dx_display,*fp_x_grid_dat;
	FILE *fp_dy,*fp_dy_display,*fp_y_grid_dat;
	FILE *fp_dz,*fp_dz_display,*fp_z_grid_dat;

	fp_grid_input = fopen(in_file,"r");
	fscanf(fp_grid_input,"%lf %lf %lf %lf",&x1,&x2,&x3,&x4);
	fscanf(fp_grid_input,"%lf %lf %lf",&dx_min,&dx_max,&ax1);
	fscanf(fp_grid_input,"%lf %lf %lf %lf",&y_1,&y2,&y3,&y4);
	fscanf(fp_grid_input,"%lf %lf %lf",&dy_min,&dy_max,&ay1);
	fscanf(fp_grid_input,"%lf %lf %lf %lf",&z1,&z2,&z3,&z4);
	fscanf(fp_grid_input,"%lf %lf %lf",&dz_min,&dz_max,&az1);
	fclose(fp_grid_input);
	
	fprintf(stderr,"%g\t%g\t%g\t%g\n",x1,x2,x3,x4);
	fprintf(stderr,"%g\t%g\t%g\n",dx_min,dx_max,ax1);
	fprintf(stderr,"%g\t%g\t%g\t%g\n",y_1,y2,y3,y4);
	fprintf(stderr,"%g\t%g\t%g\n",dy_min,dy_max,ay1);
	fprintf(stderr,"%g\t%g\t%g\t%g\n",z1,z2,z3,z4);
	fprintf(stderr,"%g\t%g\t%g\n",dz_min,dz_max,az1);
	
	sprintf(mes,"%s/x_grid_spacing",output_directory);
	fp_dx = fopen(mes,"w");
	sprintf(mes,"%s/y_grid_spacing",output_directory);
	fp_dy = fopen(mes,"w");
	sprintf(mes,"%s/z_grid_spacing",output_directory);
	fp_dz = fopen(mes,"w");
	sprintf(mes,"%s/x_grid_display",output_directory);
	fp_dx_display = fopen(mes,"w");
	sprintf(mes,"%s/y_grid_display",output_directory);
	fp_dy_display = fopen(mes,"w");
	sprintf(mes,"%s/z_grid_display",output_directory);
	fp_dz_display = fopen(mes,"w");
	sprintf(mes,"%s/x_grid.dat",output_directory);
	fp_x_grid_dat = fopen(mes,"w");
	sprintf(mes,"%s/y_grid.dat",output_directory);
	fp_y_grid_dat = fopen(mes,"w");
	sprintf(mes,"%s/z_grid.dat",output_directory);
	fp_z_grid_dat = fopen(mes,"w");

	sprintf(mes,"%s/lip_sites.dat",output_directory);
	fp_lip_info = fopen(mes,"w");

	ax2 = (dx_max / dx_min) - 1.0;
	ay2 = (dy_max / dy_min) - 1.0;
	az2 = (dz_max / dz_min) - 1.0;

	i_min = i = 0;
        x_grid[0] = x1;
        dx[0] = dx_min * stretch_factor(x_grid[0] - x2,ax1,ax2);
        while(1) {
                x_grid[i+1] = x_grid[i] + dx[i];
                if(x_grid[i+1] < x2) {
                        dx[i+1] = dx_min * stretch_factor(x_grid[i+1] - x2,ax1,ax2);
                }
                else if(x_grid[i+1] < x3) {
                        dx[i+1] = dx_min;
                }
                else if(x_grid[i+1] < x4) {
                        dx[i+1] = dx_min * stretch_factor(x_grid[i+1] - x3,ax1,ax2);
                }
                else if(x_grid[i+1] >= x4) {
                        i_max = i+1;
                        x_max = x_grid[i+1];
                        dx[i+1] = dx[i];        // not really needed
                        break;
                }
                ++i;
        }

	for(i = i_min; i <= i_max; i++) {
		fprintf(fp_dx,"%g\t%g\n",x_grid[i],dx[i]);
		fprintf(fp_dx_display,"%g\t%g\n",x_grid[i],y2);
		fprintf(fp_dx_display,"%g\t%g\n",x_grid[i],y3);
		fprintf(fp_x_grid_dat,"%d\t%g\t%g\n",i,x_grid[i],dx[i]);
	}

	j_min = j = 0;
	y_grid[0] = y_1;
	dy[0] = dy_min * stretch_factor(y_grid[0] - y2,ay1,ay2);
	while(1) {
		y_grid[j+1] = y_grid[j] + dy[j];
		if(y_grid[j+1] < y2) {
			dy[j+1] = dy_min * stretch_factor(y_grid[j+1] - y2,ay1,ay2);
		}
		else if(y_grid[j+1] < y3) {
			dy[j+1] = dy_min;
		}
		else if(y_grid[j+1] < y4) {
			dy[j+1] = dy_min * stretch_factor(y_grid[j+1] - y3,ay1,ay2);
		}
		else if(y_grid[j+1] >= y4) {
			j_max = j+1;
			y_max = y_grid[j+1];
			dy[j+1] = dy[j];	// not really needed
			break;
		}
		++j;
	}

	for(j = j_min; j <= j_max; j++) {
		fprintf(fp_dy,"%g\t%g\n",y_grid[j],dy[j]);
		fprintf(fp_dy_display,"%g\t%g\n",x2,y_grid[j]);
		fprintf(fp_dy_display,"%g\t%g\n",x3,y_grid[j]);
		fprintf(fp_y_grid_dat,"%d\t%g\t%g\n",j,y_grid[j],dy[j]);
	}

	k_min = k = 0;
	z_grid[0] = z1;
	dz[0] = dz_min * stretch_factor(z_grid[0] - z2,az1,az2);
	while(1) {
		z_grid[k+1] = z_grid[k] + dz[k];
		if(z_grid[k+1] < z2) {
			dz[k+1] = dz_min * stretch_factor(z_grid[k+1] - z2,az1,az2);
		}
		else if(z_grid[k+1] < z3) {
			dz[k+1] = dz_min;
		}
		else if(z_grid[k+1] < z4) {
			dz[k+1] = dz_min * stretch_factor(z_grid[k+1] - z3,az1,az2);
		}
		else if(z_grid[k+1] >= z4) {
			k_max = k+1;
			z_max = z_grid[k+1];
			dz[k+1] = dz[k];	// not really needed
			break;
		}
		++k;
	}

	for(k = k_min; k <= k_max; k++) {
		fprintf(fp_dz,"%g\t%g\n",z_grid[k],dz[k]);
		fprintf(fp_dz_display,"%g\t%g\n",x2,z_grid[k]);
		fprintf(fp_dz_display,"%g\t%g\n",x3,z_grid[k]);
		fprintf(fp_z_grid_dat,"%d\t%g\t%g\n",k,z_grid[k],dz[k]);
	}

	fprintf(stderr,"number of grid points = %d x %d x %d = %d\n",i_max+1,j_max+1,k_max+1,(i_max+1)*(j_max+1)*(k_max+1));
	fprintf(stderr,"grid points for uniform grid = %d x %d x %d = %d\n",(int)((x4-x1)/dx_min)+1,(int)((y4-y_1)/dy_min)+1,(int)((z4-z1)/dz_min)+1,((int)((x4-x1)/dx_min)+1)*((int)((y4-y_1)/dy_min)+1)*((int)((z4-z1)/dz_min)+1));
	fprintf(stderr,"x_max = %g\ty_max = %g\tz_max = %g\n",x_grid[i_max],y_grid[j_max],z_grid[k_max]);

	fclose(fp_dx);
	fclose(fp_dy);
	fclose(fp_dz);
	fclose(fp_dx_display);
	fclose(fp_dy_display);
	fclose(fp_dz_display);
	fclose(fp_x_grid_dat);
	fclose(fp_y_grid_dat);
	fclose(fp_z_grid_dat);

	return(0);
}

double
stretch_factor(double v, double a1, double a2)
{
	double val;
	double dum;
	double x;

	x = a1*fabs(v);
	dum = (exp(x) - exp(-x))/(exp(x) + exp(-x));
	val = 1.0 + a2 * pow(dum,2.0);

	return(val);
}		

int
place_recorder()
{
	FILE *fp_rec,*fp,*fp_u_init;
	int i,j,k,n,m,p,q,r;
        char tmp[N_Y_MAX],tmp2[N_Y_MAX];
        int i_rec,j_rec;
	int init_flag;
	char velocity_flag[100];
	int n_entries;
        char mes[800],c_junk;
	int n_layer[N_Z_MAX];
	int i_init_bare_min,i_init_bare_max;
	int j_init_bare_min,j_init_bare_max;
	int k_init_bare_min,k_init_bare_max;
	int im,jm,km;
	FILE *fp_marked;
	int i_circle_center,j_circle_center,k_circle_center;
	int r_circle;
	int i_circle,j_circle,k_circle;
	int j_offset;
	int i_start_shell,i_finish_shell;
	int circle_flag,fill_in_flag;
	int k_shell;

	int i_left_min;
	int i_left_max;
	int j_left_min;
	int j_left_max;
	int k_left_min;
	int k_left_max;

	int i_right_min;
	int i_right_max;
	int j_right_min;
	int j_right_max;
	int k_right_min;
	int k_right_max;
	int n_circles = 0;
	int j_fill_up_start,j_fill_up_end;
	int k_fill_up_start[N_Y_MAX],k_fill_up_end[N_Y_MAX];
	int k_fill_down_start[N_Y_MAX],k_fill_down_end[N_Y_MAX];

// printf("r = %d: %d %d %d\n",r_circle,i+m,j,k_circle_center+k_shell);

	int i_corner[2],j_corner[2],k_corner[2];
	double i_origin[2],j_origin[2],k_origin[2]; // origin positions of lips, LEFT or RIGHT
	int i_pivot[2],j_pivot[2],k_pivot[2];
	double i_width[2],j_width[2],k_width[2];

	sprintf(mes,"%s/marked_points",output_directory);
	fp_marked = fopen(mes,"w");

        fp = fopen(recorder_geometry,"r");
        if(fp == NULL) {
                fprintf(stderr,"can't open %s for geometry data\n",recorder_geometry);
                exit(0);
        }

        while(1) {
                if(fgets(tmp,300,fp) == NULL) {
                        fprintf(stderr, "error 0 in reading geometry file\n");
                        (void)exit(0);
                }
                if(tmp[0] == '#') break;
        }

        fgets(tmp,300,fp);
        sscanf(tmp,"%lf",&dx_bare);
        fgets(tmp,300,fp);
        sscanf(tmp,"%d",&i_max_bare);
        fgets(tmp,300,fp);
        sscanf(tmp,"%d",&j_max_bare);
        fgets(tmp,300,fp);
        sscanf(tmp,"%d",&k_max_bare);
        fgets(tmp,300,fp);
        sscanf(tmp,"%d",&i_recorder_bare);
        fgets(tmp,300,fp);
        sscanf(tmp,"%d",&j_recorder_bare);
        fgets(tmp,300,fp);
        sscanf(tmp,"%d",&k_recorder_bare);

	x_recorder_bare = i_recorder_bare * dx_bare;
	y_recorder_bare = j_recorder_bare * dx_bare;
	z_recorder_bare = k_recorder_bare * dx_bare;

	bare_geometry = allocate_array_of_chars(i_max_bare+2,j_max_bare+2,k_max_bare+2);

        for(i = 0; i <= i_max_bare; i++) {
                for(j = 0; j <= j_max_bare; j++) {
                	for(k = 0; k <= k_max_bare; k++) {
                        	bare_geometry[i][j][k] = OPEN;
			}
                }
        }
        for(i = 0; i <= i_max_bare; i++) {
                for(j = 0; j <= j_max_bare; j++) {
                	bare_geometry[i][j][0] = LOWER_Z_BOUNDARY;
                	bare_geometry[i][j][k_max_bare] = UPPER_Z_BOUNDARY;
		}
        }
        for(j = 0; j <= j_max_bare; j++) {
                for(k = 0; k <= k_max_bare; k++) {
                	bare_geometry[0][j][k] = LOWER_X_BOUNDARY;
                	bare_geometry[i_max_bare][j][k] = UPPER_X_BOUNDARY;
		}
        }
        for(i = 0; i <= i_max_bare; i++) {
                for(k = 0; k <= k_max_bare; k++) {
                	bare_geometry[i][0][k] = LOWER_Y_BOUNDARY;
                	bare_geometry[i][j_max_bare][k] = UPPER_Y_BOUNDARY;
		}
        }

        i = i_recorder_bare;
        j = j_recorder_bare;
        k = k_recorder_bare;	// k is the index of the next layer to be populated

// now read in the first "page" of the instrument from the rec_geom file
//	the pages are indexed using q starting with q = 1
// the velocity scale factor for each page will be stored in v_scale_factor_bare[q]

// each of these pages will be repeated a certain number of times creating
// the instrument layer by layer - each layer is parallel to the x-y plane

	q = 0;
        fgets(tmp,N_Y_MAX,fp);

	while(1) {
		if(tmp[0] != '@') {
            		fprintf(stderr,"error 1 in geometry file\n");
               		(void)exit(0);
        	}
		++q;
		sscanf(tmp+1,"%d %lf",&(n_layer[q]),&(v_scale_factor_bare[k]));	// n_layer[q] is the number of times to repeat this page
										// v_scale_factor_bare[k] is the scale factor for all of these layers
		if(n_layer[q] < 1) break;	// done with pages

fprintf(stderr,"page %d\tplaced starting in layer %d\trepeat %d times\tv_factor = %g\n",q,k,n_layer[q],v_scale_factor_bare[k]);
// starting a new page
// place the first layer of this page in the plane k which is initially set to k_recorder_bare
        	i = i_recorder_bare;		// each page starts at this value of i
        	while(1) {	// this loop reads in the page one line (value of i) at a time and fills out the
				// appropriate region in the k direction
               		if(fgets(tmp,N_Y_MAX,fp) == NULL) {
               	       		fprintf(stderr,"error 2 in geometry file\n");
               	        	(void)exit(0);
               	 	}
               	 	if(tmp[0] == '@') break;	// this means we are done with this page - start this loop over again for the next page

                	sscanf(tmp,"%d %s",&n,tmp2);	// n is the number of times to repeat this row (parallel to j (y) direction)
							// tmp2 contains the data for this row for which points are solid and open
fprintf(stderr,"%s\n",tmp2);
fprintf(stderr,"now repeat this row %d times starting at i = %d\n",n,i);
							// next determine if this is a regular (flat) row or if this has circular symmetry
			circle_flag = NO;
                        for(p = 0; p < strlen(tmp2); p++) {	
				if(tmp2[p] == 'C') {
					circle_flag = YES;
					break;
				}
			}

			if(circle_flag == NO) {		// this is a flat row - not circular
                		for(m = 0; m < n; m++) {
                       			j = j_recorder_bare;
                        		for(p = 0; p < strlen(tmp2); p++) {	// read the data for this row one value of j (y direction) at a time
                               		 	if(tmp2[p] == 'X') bare_geometry[i+m][j][k] = SOLID;
                               		 	if(tmp2[p] == '.') bare_geometry[i+m][j][k] = OPEN;
                               		 	if(tmp2[p] == 'I') {
							bare_geometry[i+m][j][k] = INIT_POINT;
//fprintf(stderr,"init point at %d\t%d\t%d\tv_scale = %g\n",i+m,j,k,v_scale_factor_bare[k]);
//fprintf(fp_diag,"%d\t%d\t%d\t%d\n",i,j,k,geometry[i][j][k]);
						}
                               		 	if(tmp2[p] == 'M') {	// these are special marking points
							bare_geometry[i+m][j][k] = MARKED;
						}
                               		 	if(tmp2[p] == 'L') bare_geometry[i+m][j][k] = LEFT_LIP;
                               		 	if(tmp2[p] == 'R') bare_geometry[i+m][j][k] = RIGHT_LIP;
						++j;
					}
						// now fill out the corresponding regions along the k direction
					for(p = 1; p < n_layer[q]; p++) {
						for(j = 0; j <= j_max_bare; j++) {
							bare_geometry[i+m][j][k+p] = bare_geometry[i+m][j][k];
						}
						v_scale_factor_bare[k + p] = v_scale_factor_bare[k];
					}
				}
				i += n;
			}
			else {			// deal instead with circular regions (shells)
                        	for(p = 0; p < strlen(tmp2); p++) {	
					if(tmp2[p] == 'O') {
						j_circle_center = p + j_recorder_bare;
						break;
					}
				}
fprintf(stderr,"\n\nstarting next bunch of circular regions centered at i/j = %d/%d\n",i+m,j_circle_center);
				for(m = 0; m < n; m++) {
// fprintf(stderr,"\nstarting new circular region centered at i/j = %d/%d\n",i+m,j_circle_center);
// fprintf(stderr,"%s\n",tmp2);

					fill_in_flag = FIRST_CIRCLE;	// this flag is FIRST, MIDDLE, or LAST_CIRCLE in the shell
						// now do each circle that makes up this shell
                        		for(p = 0; p < strlen(tmp2); p++) {	
						k_circle_center = k + n_layer[q]/2;
						if(tmp2[p] == 'O') break;
						if(tmp2[p] == 'C') {
							r_circle = j_circle_center - j_recorder_bare - p; // find radius of shell
							if(r_circle < 0) {
								fprintf(stderr,"error in geometry file - r_shell negative\n");
								exit(0);
							}
//							k_circle_center = k + n_layer[q]/2;
// fprintf(stderr,"fill in circular shell of radius r = %d centered on i/j/k = %d/%d/%d\n",r_circle,i+m,j_circle_center,k_circle_center);
							for(j = j_circle_center - r_circle; j <= j_circle_center + r_circle; j++) {
								k_shell = sqrt(r_circle*r_circle - (j_circle_center-j)*(j_circle_center-j));
// fprintf(stdout,"\tfilling in circular shell i/j/k %d/%d/%d\tand %d\n",i+m,j,k_circle_center+k_shell,k_circle_center-k_shell);
								bare_geometry[i+m][j][k_circle_center+k_shell] = SOLID;
								bare_geometry[i+m][j][k_circle_center-k_shell] = SOLID;

								if(fill_in_flag == FIRST_CIRCLE) {
									k_fill_up_start[j] = k_fill_up_end[j] = k_circle_center+k_shell;
									k_fill_down_start[j] = k_fill_down_end[j] = k_circle_center-k_shell;
								}
								if(fill_in_flag == LAST_CIRCLE) {
									k_fill_up_end[j] = k_circle_center+k_shell;
									k_fill_down_end[j] = k_circle_center-k_shell;
								}
// printf("r = %d: %d %d %d %d\n",r_circle,i+m,j,k_circle_center,k_shell);
							}
							if(fill_in_flag == FIRST_CIRCLE) {
								j_fill_up_start = j_circle_center - r_circle;
								j_fill_up_end = j_circle_center + r_circle;
								fill_in_flag = LAST_CIRCLE;
							}
						}
					}
// printf("\tj_fill_up start/end: %d %d\n",j_fill_up_start,j_fill_up_end);
// printf("\tj / k_fill_up_start / end %d %d %d\n",j_fill_up_start,k_fill_up_start[j_fill_up_start],k_fill_up_end[j_fill_up_start]);
// printf("\tj / k_fill_up_start / end %d %d %d\n",(j_fill_up_end+j_fill_up_start)/2,k_fill_up_start[(j_fill_up_end+j_fill_up_start)/2],k_fill_up_end[(j_fill_up_end+j_fill_up_start)/2]);

					for(j = j_fill_up_start; j <= j_fill_up_end; j++) {
						for(r = k_fill_up_start[j]; r >= k_fill_up_end[j]; r--) {
							bare_geometry[i+m][j][r] = SOLID;
// printf("\tfill: %d %d %d\n",i+m,j,r);
						}
// printf("\n");
						for(r = k_fill_down_start[j]; r <= k_fill_down_end[j]; r++) {
							bare_geometry[i+m][j][r] = SOLID;
// printf("\tfill: %d %d %d\n",i+m,j,r);
						}
// printf("\n");
					}
// exit(0);
				}
				i += n;
			}
		}
		k += n_layer[q];
	}

fprintf(stderr,"done with inserting pages in bare_geometry\n");

// now find the initialization points in the channel
	i_init_bare_min = i_max_bare;
	i_init_bare_max = 0;
	j_init_bare_min = j_max_bare;
	j_init_bare_max = 0;
	k_init_bare_min = k_max_bare;
	k_init_bare_max = 0;
fprintf(fp_diag,"init points in bare and final coords\n");
        for(i = 0; i <= i_max_bare; i++) {
                for(j = 0; j <= j_max_bare; j++) {
                	for(k = 0; k <= k_max_bare; k++) {
                        	if(bare_geometry[i][j][k] == INIT_POINT) {
					if(i < i_init_bare_min) i_init_bare_min = i;
					if(i > i_init_bare_max) i_init_bare_max = i;
					if(j < j_init_bare_min) j_init_bare_min = j;
					if(j > j_init_bare_max) j_init_bare_max = j;
					if(k < k_init_bare_min) k_init_bare_min = k;
					if(k > k_init_bare_max) k_init_bare_max = k;
fprintf(fp_diag,"%d\t%d\t%d\t%d\t%d\t%d\tv_scale = %g\n",i,i_trans(i),j,j_trans(j),k,k_trans(k),v_scale_factor_bare[k]);
// fprintf(fp_diag,"k_recorder_bare = %d\n",k_recorder_bare);
fflush(fp_diag);
// exit(0);
				}
			}
                }
        }
fprintf(fp_diag,"done with init points in bare coords\n");
fprintf(stderr,"done with init points in bare coords\n");

	i_init_min = i_trans(i_init_bare_min);
	i_init_max = i_trans(i_init_bare_max);
	j_init_min = j_trans(j_init_bare_min);
	j_init_max = j_trans(j_init_bare_max);
	k_init_min = k_trans(k_init_bare_min);
	k_init_max = k_trans(k_init_bare_max);

fprintf(fp_diag,"init points in channel\n");
fprintf(fp_diag,"%d-%d\t%d-%d\t%d-%d\n",
i_init_bare_min,i_init_bare_max,j_init_bare_min,j_init_bare_max,k_init_bare_min,k_init_bare_max);

// the next line contains info for u_init.dat
	fgets(tmp,600,fp);
	sprintf(tmp2,"%s/u_init.dat",output_directory);
	fp_u_init = fopen(tmp2,"w");
	fprintf(fp_u_init,"%s\n",tmp);
	fclose(fp_u_init);

// places where sound is recorded
//        fscanf(fp,"%d %d %d",&i_record_1_bare,&j_record_1_bare,&k_record_1_bare);
//        fscanf(fp,"%d %d %d",&i_record_2_bare,&j_record_2_bare,&k_record_2_bare);
//        fscanf(fp,"%d %d %d",&i_record_3_bare,&j_record_3_bare,&k_record_3_bare);

	i = 0;
	while(1) {
		velocity_flag[0] = '\0';
		++i;
		if(fgets(tmp,300,fp) == NULL) break;
fprintf(stderr,"%s\n",tmp);
        	n_entries = sscanf(tmp,"%d %d %d %s",&(i_record_bare[i]),&(j_record_bare[i]),&(k_record_bare[i]),velocity_flag);
fprintf(stderr,"%d %d %d %s\n",i_record_bare[i],j_record_bare[i],k_record_bare[i],velocity_flag);
		n_record_flag[i] = SOUND;
		if((n_entries == 4) && (velocity_flag[0] == 'v')) {
			n_record_flag[i] = VELOCITY;
			fprintf(stderr,"record velocity at location %d (%d / %d / %d)\n",i,i_record_bare[i],j_record_bare[i],k_record_bare[i]);
		}
	}
	n_record_points = i - 1;

        fclose(fp);

//	next extract lip information from bare_geometry[][][] and put it into file with fp_lips
fprintf(stderr,"start extracting lip info\n");

	i_left_min = i_max_bare; 
	i_left_max = 0;
	j_left_min = j_max_bare;
	j_left_max = 0;
	k_left_min = k_max_bare;
	k_left_max = 0;

	i_right_min = i_max_bare;
	i_right_max = 0;
	j_right_min = j_max_bare;
	j_right_max = 0;
	k_right_min = k_max_bare;
	k_right_max = 0;

        for(i = 0; i <= i_max_bare; i++) {
                for(j = 0; j <= j_max_bare; j++) {
                	for(k = 0; k <= k_max_bare; k++) {
                        	if(bare_geometry[i][j][k] == LEFT_LIP) {
					if(i < i_left_min) i_left_min = i;
					if(i > i_left_max) i_left_max = i;
					if(j < j_left_min) j_left_min = j;
					if(j > j_left_max) j_left_max = j;
					if(k < k_left_min) k_left_min = k;
					if(k > k_left_max) k_left_max = k;
				}
                        	if(bare_geometry[i][j][k] == RIGHT_LIP) {
					if(i < i_right_min) i_right_min = i;
					if(i > i_right_max) i_right_max = i;
					if(j < j_right_min) j_right_min = j;
					if(j > j_right_max) j_right_max = j;
					if(k < k_right_min) k_right_min = k;
					if(k > k_right_max) k_right_max = k;
				}
			}
		}
	}

fprintf(stderr,"lip boundaries: left: %d/%d\t%d/%d\t%d/%d\n",i_left_min,i_left_max,j_left_min,j_left_max,k_left_min,k_left_max);
fprintf(stderr,"lip boundaries: right: %d/%d\t%d/%d\t%d/%d\n",i_right_min,i_right_max,j_right_min,j_right_max,k_right_min,k_right_max);

	i_corner[LEFT] = i_left_min;
	j_corner[LEFT] = j_left_min;
	k_corner[LEFT] = k_left_min;
	i_corner[RIGHT] = i_right_min;
	j_corner[RIGHT] = j_right_min;
	k_corner[RIGHT] = k_right_min;

	i_pivot[LEFT] = i_left_min;
	j_pivot[LEFT] = j_left_min - 1;

	i_pivot[RIGHT] = i_right_min;
	j_pivot[RIGHT] = j_right_max + 1;

	i_origin[LEFT] = ((double)(i_left_max + i_left_min)) / 2;
	j_origin[LEFT] = ((double)(j_left_max + j_left_min)) / 2;

	i_origin[RIGHT] = ((double)(i_right_max + i_right_min)) / 2;
	j_origin[RIGHT] = ((double)(j_right_max + j_right_min)) / 2;

	i_width[LEFT] = i_left_max - i_origin[LEFT];
	j_width[LEFT] = j_left_max - j_origin[LEFT];

	i_width[RIGHT] = i_right_max - i_origin[RIGHT];
	j_width[RIGHT] = j_right_max - j_origin[RIGHT];

fprintf(stderr,"lip info: left: %d/%d\t%g/%g\t%g/%g\n",i_pivot[LEFT],j_pivot[LEFT],i_origin[LEFT],j_origin[LEFT],i_width[LEFT],j_width[LEFT]);
fprintf(stderr,"lip info: right: %d/%d\t%g/%g\t%g/%g\n",i_pivot[RIGHT],j_pivot[RIGHT],i_origin[RIGHT],j_origin[RIGHT],i_width[RIGHT],j_width[RIGHT]);

        for(i = i_left_min - 2; i <= i_left_max + 2; i++) {
        	for(j = j_left_min - 2; j <= j_right_max + 2; j++) {
                	k = (k_left_max + k_left_min) / 2;
			if((i == i_pivot[LEFT]) && (j == j_pivot[LEFT])) {
				fprintf(fp_lip_info,"P");
			}
			else if((i == i_pivot[RIGHT]) && (j == j_pivot[RIGHT])) {
				fprintf(fp_lip_info,"P");
			}
			else if((i == (int)i_origin[LEFT]) && (j == (int)j_origin[LEFT])) {
				fprintf(fp_lip_info,"O");
			}
			else if((i == (int)i_origin[RIGHT]) && (j == (int)j_origin[RIGHT])) {
				fprintf(fp_lip_info,"O");
			}
			else {
                        	if(bare_geometry[i][j][k] == LEFT_LIP) fprintf(fp_lip_info,"L");
                        	if(bare_geometry[i][j][k] == RIGHT_LIP) fprintf(fp_lip_info,"R");
                        	if(bare_geometry[i][j][k] == OPEN) fprintf(fp_lip_info,".");
                        	if(bare_geometry[i][j][k] == SOLID) fprintf(fp_lip_info,"X");
			}
		}
		fprintf(fp_lip_info,"\n");
	}
fflush(fp_lip_info);
fprintf(stderr,"done with lip boundaries\n");

//	now convert from the bare geometry (uniform grid) to the new geometry (nonuniform grid)
//	will have to adjust the placement of some points but keep them close to the original
//	locations in terms of real space (x and y)

//	first locate the grid point closest to the desired location of the lower left
//	point on the recorder

	for(i = 0; i < i_max; i++) {
		if(fabs(x_grid[i] - x_recorder_bare) <= dx[i]/2.0) {
			i_recorder_corner = i;
			break;
		}
	}
	for(j = 0; j < j_max; j++) {
		if(fabs(y_grid[j] - y_recorder_bare) <= dy[j]/2.0) {
			j_recorder_corner = j;
			break;
		}
	}
	for(k = 0; k < k_max; k++) {
		if(fabs(z_grid[k] - z_recorder_bare) <= dz[k]/2.0) {
			k_recorder_corner = k;
			break;
		}
	}

//	now transfer the recorder location in bare_geometry[][][] to points in geometry[][][]
//	have to be careful not to overwrite points that should be open with solid layers
//	that "leak" over
//	new code on 5-22-16 to deal with this using the TMP_OPEN state

	geometry = allocate_array_of_chars(i_max+2,j_max+2,k_max+2);

	for(i = 0; i < i_max; i++) {
		for(j = 0; j < j_max; j++) {
			for(k = 0; k < k_max; k++) {
				geometry[i][j][k] = TMP_OPEN;
			}
		}
	}
	for(i = 0; i <= i_max; i++) {
		for(k = 0; k < k_max; k++) {
                	geometry[i][0][k] = LOWER_Y_BOUNDARY;
                	geometry[i][j_max][k] = UPPER_Y_BOUNDARY;
		}
        }
        for(j = 0; j <= j_max; j++) {
		for(k = 0; k < k_max; k++) {
                	geometry[0][j][k] = LOWER_X_BOUNDARY;
                	geometry[i_max][j][k] = UPPER_X_BOUNDARY;
		}
        }
        for(i = 0; i <= i_max; i++) {
        	for(j = 0; j <= j_max; j++) {
                		geometry[i][j][0] = LOWER_Z_BOUNDARY;
                		geometry[i][j][k_max] = UPPER_Z_BOUNDARY;
		}
        }

fprintf(fp_diag,"converting\n");
// this is the new version - this should get rid of over-writing problem noted in previous versions
// deal with marked points first
	for(i = 0; i < i_max_bare; i++) {
		for(j = 0; j < j_max_bare; j++) {
			for(k = 0; k < k_max_bare; k++) {
				if(bare_geometry[i][j][k] == MARKED) {
					bare_geometry[i][j][k] = OPEN;
					im = i_trans(i);
					jm = j_trans(j);
					km = k_trans(k);
					fprintf(fp_marked,"%d\t%d\t%d\t%g\t%g\t%g\n",im,jm,km,x_grid[im],y_grid[jm],z_grid[km]);
// fprintf(stderr,"found a marked point\n");
// exit(0);
				}
			}
		}
	}
	fclose(fp_marked);

	for(i = 0; i < i_max_bare; i++) {
		for(j = 0; j < j_max_bare; j++) {
			for(k = 0; k < k_max_bare; k++) {
				if((bare_geometry[i][j][k] == SOLID) && (geometry[i_trans(i)][j_trans(j)][k_trans(k)] == TMP_OPEN)) {
					geometry[i_trans(i)][j_trans(j)][k_trans(k)] = SOLID;
				}
				if(bare_geometry[i][j][k] == OPEN) {
					geometry[i_trans(i)][j_trans(j)][k_trans(k)] = OPEN;
				}
				if(bare_geometry[i][j][k] == LEFT_LIP) {
					geometry[i_trans(i)][j_trans(j)][k_trans(k)] = LEFT_LIP;
				}
				if(bare_geometry[i][j][k] == RIGHT_LIP) {
					geometry[i_trans(i)][j_trans(j)][k_trans(k)] = RIGHT_LIP;
				}
			}
		}
	}
	for(i = 0; i < i_max_bare; i++) {
		for(j = 0; j < j_max_bare; j++) {
			for(k = 0; k < k_max_bare; k++) {
				if(bare_geometry[i][j][k] == INIT_POINT) {
					geometry[i_trans(i)][j_trans(j)][k_trans(k)] = INIT_POINT;
					v_scale_factor[k_trans(k)] = v_scale_factor_bare[k];
fprintf(fp_diag,"%d\t%d\t%d\t%d\t%d\t%d\n",i,i_trans(i),j,j_trans(j),k,k_trans(k));
				}
			}
		}
	}
fprintf(fp_diag,"finished converting\n");

fprintf(fp_diag,"checking\n");
for(i = i_init_min; i <= i_init_max; i++) {
	for(j = j_init_min; j <= j_init_max; j++) {
		for(k = k_init_min; k <= k_init_max; k++) {
			if(geometry[i][j][k] != INIT_POINT) {
fprintf(fp_diag,"%d\t%d\t%d\n",i,j,k);
			}
		}
	}
}
fprintf(fp_diag,"done checking\n");

//        i_record_1 = i_trans(i_record_1_bare);
//	j_record_1 = j_trans(j_record_1_bare);
//	k_record_1 = k_trans(k_record_1_bare);
//       i_record_2 = i_trans(i_record_2_bare);
//	j_record_2 = j_trans(j_record_2_bare);
//	k_record_2 = k_trans(k_record_2_bare);
//        i_record_3 = i_trans(i_record_3_bare);
//	j_record_3 = j_trans(j_record_3_bare);
//	k_record_3 = k_trans(k_record_3_bare);

	for(i = 1; i <= n_record_points; i++) {
        	i_record[i] = i_trans(i_record_bare[i]);
		j_record[i] = j_trans(j_record_bare[i]);
		k_record[i] = k_trans(k_record_bare[i]);
	}
        
fprintf(fp_lip_info,"\n\nnow in transformed coordinates:\n");
        for(i = i_trans(i_left_min) - 2; i <= i_trans(i_left_max) + 2; i++) {
        	for(j = j_trans(j_left_min) - 2; j <= j_trans(j_right_max) + 2; j++) {
                	k = k_trans((k_left_max + k_left_min) / 2);
			if((i == i_trans(i_pivot[LEFT])) && (j == j_trans(j_pivot[LEFT]))) {
				fprintf(fp_lip_info,"P");
			}
			else if((i == i_trans(i_pivot[RIGHT])) && (j == j_trans(j_pivot[RIGHT]))) {
				fprintf(fp_lip_info,"P");
			}
			else if((i == i_trans((int)i_origin[LEFT])) && (j == j_trans((int)j_origin[LEFT]))) {
				fprintf(fp_lip_info,"O");
			}
			else if((i == i_trans((int)i_origin[RIGHT])) && (j == j_trans((int)j_origin[RIGHT]))) {
				fprintf(fp_lip_info,"O");
			}
			else {
                        	if(geometry[i][j][k] == LEFT_LIP) fprintf(fp_lip_info,"L");
                        	if(geometry[i][j][k] == RIGHT_LIP) fprintf(fp_lip_info,"R");
                        	if(geometry[i][j][k] == OPEN) fprintf(fp_lip_info,".");
                        	if(geometry[i][j][k] == SOLID) fprintf(fp_lip_info,"X");
			}
		}
		fprintf(fp_lip_info,"\n");
	}

	for(i = 0; i < i_max; i++) {
		for(j = 0; j < j_max; j++) {
			for(k = 0; k < k_max; k++) {
				if(geometry[i][j][k] == TMP_OPEN) geometry[i][j][k] = OPEN;
			}
		}
	}
        return(0);
}

/*	
 *	functions to convert from bare coordinates to coordinates on the nonuniform grid
 */

i_trans(i_bare)
int i_bare;
{
	int i;
	double x_bare;
	static int i_last = 0;

	x_bare = i_bare * dx_bare;

	if(fabs(x_bare - x_grid[0]) <= dx[0]/2.0) {
		i_last = 0;
		return(i_last);
	}
	for(i = i_last; i < i_max; i++) {
		if((fabs(x_grid[i] - x_bare) <= dx[i-1]/2.0) && (i > 0)) {
			i_last = i;
			return(i_last);
		}
		else if(fabs(x_bare - x_grid[i]) <= dx[i]/2.0) {
			i_last = i;
			return(i_last);
		}
	}
	if(fabs(x_grid[i_max] - x_bare) <= dx[i_max-1]/2.0) {
		i_last = i_max;
		return(i_last);
	}
	if(fabs(x_bare - x_grid[i_max]) <= dx[i_max]/2.0) {
		i_last = i_max;
		return(i_last);
	}
	for(i = 0; i < i_last+1; i++) {
		if(fabs(x_grid[i] - x_bare) <= dx[i-1]/2.0) {
			i_last = i;
			return(i_last);
		}
		else if(fabs(x_bare - x_grid[i]) <= dx[i]/2.0) {
			i_last = i;
			return(i_last);
		}
	}
	fprintf(stderr,"can't locate the new x coordinate\n");
	fprintf(stderr,"i_bare = %d\tx_bare = %g\n",i_bare,x_bare);
	fprintf(fp_log,"can't locate the new x coordinate\n");
	fprintf(fp_log,"i_bare = %d\tx_bare = %g\n",i_bare,x_bare);
	fflush(fp_log);
	exit(0);
	
	return(i_last);
}

j_trans(j_bare)
int j_bare;
{
	int j;
	double y_bare;
	static int j_last = 0;

	y_bare = j_bare * dx_bare;	// bare coordinate of the corner of the location to be converted

	if(fabs(y_bare - y_grid[0]) <= dy[0]/2.0) {
		j_last = 0;
		return(j_last);
	}
	for(j = j_last; j < j_max; j++) {
		if((fabs(y_grid[j] - y_bare) <= dy[j-1]/2.0) && (j > 0)) {
			j_last = j;
			return(j_last);
		}
		else if(fabs(y_bare - y_grid[j]) <= dy[j]/2.0) {
			j_last = j;
			return(j_last);
		}
	}
	if(fabs(y_grid[j_max] - y_bare) <= dy[j_max-1]/2.0) {
		j_last = j_max;
		return(j_last);
	}
	if(fabs(y_bare - y_grid[j_max]) <= dy[j_max]/2.0) {
		j_last = j_max;
		return(j_last);
	}
	for(j = 1; j < j_last+1; j++) {
		if(fabs(y_grid[j] - y_bare) <= dy[j-1]/2.0) {
			j_last = j;
			return(j_last);
		}
		else if(fabs(y_bare - y_grid[j]) <= dy[j]/2.0) {
			j_last = j;
			return(j_last);
		}
	}
	fprintf(stderr,"can't locate the new y coordinate\n");
	fprintf(stderr,"j_bare = %d\ty_bare = %g\n",j_bare,y_bare);
	fprintf(fp_log,"can't locate the new y coordinate\n");
	fprintf(fp_log,"j_bare = %d\ty_bare = %g\n",j_bare,y_bare);
	fflush(fp_log);
	exit(0);
	
	return(0);
}

k_trans(k_bare)
int k_bare;
{
	int k;
	double z_bare;
	static int k_last = 0;

	z_bare = k_bare * dx_bare;

// fprintf(fp_log,"k_bare = %d\tz_bare = %g\tk_last = %d\n",k_bare,z_bare,k_last);
// fprintf(fp_log,"z_grid[0] = %g\tdz[0] = %g\n",z_grid[0],dz[0]);

	if(fabs(z_bare - z_grid[0]) <= dz[0]/2.0) {
		k_last = 0;
		return(k_last);
	}
	for(k = k_last; k < k_max; k++) {
		if((fabs(z_grid[k] - z_bare) <= dz[k-1]/2.0) && (k > 0)) {
			k_last = k;
// fprintf(fp_log,"k_trans 1: k = %d\tz_grid[k] = %g\tdz[k-1] = %g\n",k,z_grid[k],dz[k-1]);
			return(k_last);
		}
		else if(fabs(z_bare - z_grid[k]) <= dz[k]/2.0) {
			k_last = k;
// fprintf(fp_log,"k_trans 2: k = %d\tz_grid[k] = %g\tdz[k] = %g\n",k,z_grid[k],dz[k]);
			return(k_last);
		}
	}
	if(fabs(z_grid[k_max] - z_bare) <= dz[k_max-1]/2.0) {
		k_last = k_max;
		return(k_last);
	}
	if(fabs(z_bare - z_grid[k_max]) <= dz[k_max]/2.0) {
		k_last = k_max;
		return(k_last);
	}
	for(k = 1; k <= k_last; k++) {
		if(fabs(z_grid[k] - z_bare) <= dz[k-1]/2.0) {
			k_last = k;
			return(k_last);
		}
		else if(fabs(z_bare - z_grid[k]) <= dz[k]/2.0) {
			k_last = k;
			return(k_last);
		}
	}
	fprintf(stderr,"can't locate the new z coordinate\n");
	fprintf(stderr,"k_bare = %d\tz_bare = %g\n",k_bare,z_bare);
	fprintf(fp_log,"can't locate the new z coordinate\n");
	fprintf(fp_log,"k_bare = %d\tz_bare = %g\n",k_bare,z_bare);
	fflush(fp_log);
	exit(0);
	
	return(0);
}

/*
 *	store geometry data in a vtk file for view with paraview
 */
make_vtk_file1()
{
	int i,j,k;
	FILE *fp_vtk;
	long int n_points;
	int n_rec_points;
	int i_rec_min,j_rec_min,k_rec_min;
	int i_rec_max,j_rec_max,k_rec_max;
	char mes[200];

	i_rec_min = i_max;
	j_rec_min = j_max;
	k_rec_min = k_max;
	i_rec_max = j_rec_max = k_rec_max = 0;
	for(i = 0; i <= i_max; i++) {
		for(j = 0; j <= j_max; j++) {
			for(k = 0; k <= k_max; k++) {
				if((geometry[i][j][k] == SOLID) 
					|| (geometry[i][j][k] == INIT_POINT)) {
					if(i < i_rec_min) i_rec_min = i;
					if(i > i_rec_max) i_rec_max = i;
					if(j < j_rec_min) j_rec_min = j;
					if(j > j_rec_max) j_rec_max = j;
					if(k < k_rec_min) k_rec_min = k;
					if(k > k_rec_max) k_rec_max = k;
				}
			}
		}
	}

	n_points = (i_rec_max-i_rec_min+1)*(j_rec_max-j_rec_min+1)*(k_rec_max-k_rec_min+1);

	sprintf(mes,"%s/geometry.vtk",vtk_directory);
	fp_vtk = fopen(mes,"w");

	fprintf(fp_vtk,"# vtk DataFile Version 2.0\n");
	fprintf(fp_vtk,"Recorder Dataset\n");
	fprintf(fp_vtk,"ASCII\n");
	fprintf(fp_vtk,"DATASET STRUCTURED_GRID\n");
	fprintf(fp_vtk,"DIMENSIONS %d %d %d\n",k_rec_max-k_rec_min+1,j_rec_max-j_rec_min+1
									,i_rec_max-i_rec_min+1);
									
	fprintf(fp_vtk,"POINTS %ld float\n",n_points);

	for(i = i_rec_min; i <= i_rec_max; i++) {
		for(j = j_rec_min; j <= j_rec_max; j++) {
			for(k = k_rec_min; k <= k_rec_max; k++) {
//				fprintf(fp_vtk,"%g\t%g\t%g\n",dx_min*i,dy_min*j,dz_min*k);
				fprintf(fp_vtk,"%g\t%g\t%g\n",x_grid[i],y_grid[j],z_grid[k]);
			}
		}
	}

	fprintf(fp_vtk,"\n");
	fprintf(fp_vtk,"POINT_DATA %ld\n",n_points);
	fprintf(fp_vtk,"SCALARS rho float 1\n");
	fprintf(fp_vtk,"LOOKUP_TABLE default\n");

	for(i = i_rec_min; i <= i_rec_max; i++) {
		for(j = j_rec_min; j <= j_rec_max; j++) {
			for(k = k_rec_min; k <= k_rec_max; k++) {
				if(geometry[i][j][k] == SOLID) {
					fprintf(fp_vtk,"10.0\n");
				}
				else if(geometry[i][j][k] == INIT_POINT) {
					fprintf(fp_vtk,"5.0\n");
				}
//				else if((i == i_record_1) && (j == j_record_1) && (k == k_record_1)) {
//					fprintf(fp_vtk,"0.25\n");
//				}
//				else if((i == i_record_2) && (j == j_record_2) && (k == k_record_2)) {
//					fprintf(fp_vtk,"0.25\n");
//				}
//				else if((i == i_record_3) && (j == j_record_3) && (k == k_record_3)) {
//					fprintf(fp_vtk,"0.25\n");
//				}
				else {
					fprintf(fp_vtk,"0.0\n");
				}
			}
		}
	}

	fclose(fp_vtk);

	return(0);
}


/*
 *      save the geometry of the simulation:
 *	the solid portion of the recorder
 *      the points where the velocity is initialized (and held fixed) 
 *      the points where rho-vs-time is recordeded 
 */
int
save_geometry()
{
        int i,j,k;
        FILE *fp_solid,*fp_initialize,*fp_sound,*fp_instrument,*fp_velocity;
	FILE *fp_left_lip,*fp_right_lip;
        char mes[400],tmp[200];

	sprintf(mes,"%s/geometry_solid.dat",output_directory);
        fp_solid = fopen(mes,"w");
	sprintf(mes,"%s/geometry_initialize.dat",output_directory);
        fp_initialize = fopen(mes,"w");
	sprintf(mes,"%s/geometry_sound.dat",output_directory);
        fp_sound = fopen(mes,"w");
	sprintf(mes,"%s/geometry_velocity.dat",output_directory);
        fp_velocity = fopen(mes,"w");
	sprintf(mes,"%s/geometry_instrument.dat",output_directory);
        fp_instrument = fopen(mes,"w");
	sprintf(mes,"%s/geometry_left_lip.dat",output_directory);
        fp_left_lip = fopen(mes,"w");
	sprintf(mes,"%s/geometry_right_lip.dat",output_directory);
        fp_right_lip = fopen(mes,"w");

fprintf(fp_diag,"saving geom\n");
        for(i = 0; i <= i_max; i++) {
                for(j = 0; j <= j_max; j++) {
                	for(k = 0; k <= k_max; k++) {
                       		if(geometry[i][j][k] == INIT_POINT) {
                               		fprintf(fp_initialize,"%d\t%d\t%d\t%g\n",i,j,k,v_scale_factor[k]);
fprintf(fp_diag,"%d\t%d\t%d\t%d\t%g\n",i,j,k,geometry[i][j][k],v_scale_factor[k]);
				}
                        	else if((geometry[i][j][k] !=  OPEN) && (geometry[i][j][k] != MARKED)) {
                                	fprintf(fp_solid,"%d\t%d\t%d\n",i,j,k);
				}
                        }
		}
	}
        for(i = 1; i < i_max; i++) {
                for(j = 1; j < j_max; j++) {
                	for(k = 1; k < k_max; k++) {
// 5-20-15                        	if(geometry[i][j][k] !=  OPEN) {
                        	if((geometry[i][j][k] !=  OPEN) && (geometry[i][j][k] !=  INIT_POINT) && (geometry[i][j][k] != MARKED)) {
                                	fprintf(fp_instrument,"%d\t%d\t%d\n",i,j,k);
				}
                        }
		}
	}
        for(i = 1; i < i_max; i++) {
                for(j = 1; j < j_max; j++) {
                	for(k = 1; k < k_max; k++) {
	                       	if(geometry[i][j][k] == LEFT_LIP) {
                               		fprintf(fp_left_lip,"%d\t%d\t%d\n",i,j,k);
				}
	                       	if(geometry[i][j][k] == RIGHT_LIP) {
                               		fprintf(fp_right_lip,"%d\t%d\t%d\n",i,j,k);
				}
                        }
		}
	}

fprintf(fp_diag,"finished saving geom\n");

	for(i = 1; i <= n_record_points; i++) {
		fprintf(fp_sound,"%d\t%d\t%d\n",i_record[i],j_record[i],k_record[i]);
fprintf(stderr,"%d\t%d\n",i,n_record_flag[i]);
		if(n_record_flag[i] == VELOCITY) fprintf(fp_velocity,"%d\t%d\t%d\n",i_record[i],j_record[i],k_record[i]);
	}

	fclose(fp_solid);
	fclose(fp_initialize);
	fclose(fp_sound);
	fclose(fp_velocity);
	fclose(fp_instrument);
	fclose(fp_left_lip);
	fclose(fp_right_lip);

	return(0);
}

/*
 *	store bare geometry data in a vtk file for view with paraview
 */
make_bare_vtk_file()
{
	int i,j,k;
	FILE *fp_vtk;
	long int n_points;
	int n_rec_points;
	int i_rec_min,j_rec_min,k_rec_min;
	int i_rec_max,j_rec_max,k_rec_max;
	char mes[200];

	i_rec_min = i_max_bare;
	j_rec_min = j_max_bare;
	k_rec_min = k_max_bare;
	i_rec_max = j_rec_max = k_rec_max = 0;
	for(i = 0; i <= i_max; i++) {
		for(j = 0; j <= j_max; j++) {
			for(k = 0; k <= k_max; k++) {
				if((bare_geometry[i][j][k] == SOLID) 
					|| (bare_geometry[i][j][k] == INIT_POINT)) {
					if(i < i_rec_min) i_rec_min = i;
					if(i > i_rec_max) i_rec_max = i;
					if(j < j_rec_min) j_rec_min = j;
					if(j > j_rec_max) j_rec_max = j;
					if(k < k_rec_min) k_rec_min = k;
					if(k > k_rec_max) k_rec_max = k;
				}
			}
		}
	}

	n_points = (i_rec_max-i_rec_min+1)*(j_rec_max-j_rec_min+1)*(k_rec_max-k_rec_min+1);

	sprintf(mes,"%s/bare_geometry.vtk",vtk_directory);
	fp_vtk = fopen(mes,"w");

	fprintf(fp_vtk,"# vtk DataFile Version 2.0\n");
	fprintf(fp_vtk,"Recorder Dataset\n");
	fprintf(fp_vtk,"ASCII\n");
	fprintf(fp_vtk,"DATASET STRUCTURED_GRID\n");
	fprintf(fp_vtk,"DIMENSIONS %d %d %d\n",k_rec_max-k_rec_min+1,j_rec_max-j_rec_min+1
									,i_rec_max-i_rec_min+1);
									
	fprintf(fp_vtk,"POINTS %ld float\n",n_points);

	for(i = i_rec_min; i <= i_rec_max; i++) {
		for(j = j_rec_min; j <= j_rec_max; j++) {
			for(k = k_rec_min; k <= k_rec_max; k++) {
				fprintf(fp_vtk,"%g\t%g\t%g\n",dx_min*i,dy_min*j,dz_min*k);
			}
		}
	}

	fprintf(fp_vtk,"\n");
	fprintf(fp_vtk,"POINT_DATA %ld\n",n_points);
	fprintf(fp_vtk,"SCALARS rho float 1\n");
	fprintf(fp_vtk,"LOOKUP_TABLE default\n");

	for(i = i_rec_min; i <= i_rec_max; i++) {
		for(j = j_rec_min; j <= j_rec_max; j++) {
			for(k = k_rec_min; k <= k_rec_max; k++) {
				if(bare_geometry[i][j][k] == SOLID) {
					fprintf(fp_vtk,"10.0\n");
				}
				else if(bare_geometry[i][j][k] == INIT_POINT) {
					fprintf(fp_vtk,"5.0\n");
				}
				else {
					fprintf(fp_vtk,"0.0\n");
				}
			}
		}
	}

	fclose(fp_vtk);

	return(0);
}

check_init_region()
{
	int i,j,k;

	fprintf(fp_diag,"checking init region\n");
	for(i = i_init_min; i <= i_init_max; i++) {
		for(j = j_init_min; j <= j_init_max; j++) {
			for(k = k_init_min; k <= k_init_max; k++) {
//				if(geometry[i][j][k] != INIT_POINT) {
					fprintf(fp_diag,"%d\t%d\t%d\t%d\n",i,j,k,geometry[i][j][k]);
//				}
			}
		}
	}
	fprintf(fp_diag,"done checking init region\n");

	return(0);
}

/*
 *	store geometry data in a second vtk file for viewing with paraview
 *	this file has only the solid points
 */
int
make_vtk_file2()
{
	int i,j,k;
	FILE *fp_vtk;
	long int n_points;
	int n_rec_points;
	int i_rec_min,j_rec_min,k_rec_min;
	int i_rec_max,j_rec_max,k_rec_max;
	char mes[200];

	i_rec_min = i_max;
	j_rec_min = j_max;
	k_rec_min = k_max;
	i_rec_max = j_rec_max = k_rec_max = 0;
	n_points = 0;
	for(i = 0; i <= i_max; i++) {
		for(j = 0; j <= j_max; j++) {
			for(k = 0; k <= k_max; k++) {
				if(geometry[i][j][k] == SOLID) {
					++n_points; 
				}
			}
		}
	}

	sprintf(mes,"%s/geometry_2.vtk",vtk_directory);
	fp_vtk = fopen(mes,"w");

	fprintf(fp_vtk,"# vtk DataFile Version 2.0\n");
	fprintf(fp_vtk,"Recorder Dataset\n");
	fprintf(fp_vtk,"ASCII\n");
	fprintf(fp_vtk,"DATASET UNSTRUCTURED_GRID\n");
//	fprintf(fp_vtk,"DIMENSIONS %d %d %d\n",k_rec_max-k_rec_min+1,j_rec_max-j_rec_min+1,i_rec_max-i_rec_min+1);
									
	fprintf(fp_vtk,"POINTS %ld float\n",n_points);

	for(i = 0; i <= i_max; i++) {
		for(j = 0; j <= j_max; j++) {
			for(k = 0; k <= k_max; k++) {
				if(geometry[i][j][k] == SOLID) 	
					 fprintf(fp_vtk,"%g\t%g\t%g\n",x_grid[i],y_grid[j],z_grid[k]);
			}
		}
	}

	fprintf(fp_vtk,"\n");
	fprintf(fp_vtk,"POINT_DATA %ld\n",n_points);
	fprintf(fp_vtk,"SCALARS rho float 1\n");
	fprintf(fp_vtk,"LOOKUP_TABLE default\n");

	for(i = 0; i <= i_max; i++) {
		for(j = 0; j <= j_max; j++) {
			for(k = 0; k <= k_max; k++) {
				if(geometry[i][j][k] == SOLID) {
					fprintf(fp_vtk,"1.0\n");
				}
			}
		}
	}

	fclose(fp_vtk);

	return(0);
}

/* 
 *  *      routine to allocate a 3d array of characters
 *   */

char ***allocate_array_of_chars(int nx, int ny, int nz)
{
        char ***a;
        int i,j;

        a = (char ***)malloc(nx * sizeof(char **));
        if(a == NULL) {
                fprintf(stderr,"first malloc failed");
                exit(0);
        }
        a[0] = (char **)malloc(nx * ny * sizeof(char *));
        if(a[0] == NULL) {
                fprintf(stderr,"second (medium size) malloc failed");
                exit(0);
        }
        a[0][0] = (char *)malloc(nx * ny * nz * sizeof(char));
        if(a[0][0] == NULL) {
                fprintf(stderr,"third (big) malloc failed");
                exit(0);
        }
        for(j = 1; j < ny; j++) {
                a[0][j] = a[0][j-1] + nz;
        }
        for(i = 1; i < nx; i++) {
                a[i] = a[i-1] + ny;
                a[i][0] = a[i-1][ny - 1] + nz;
                for(j = 1; j < ny; j++) {
                        a[i][j] = a[i][j-1] + nz;
                }
        }

	return(a);
}
